%% Calculate the log likelihood for a given nu
function log_likelihood = calculate_log_likelihood(nu, hat_theta, rep_base_se2)

% INPUT:
%   nu - Current value of nu
%   hat_theta, rep_base_se2: baseline studies estimates and SE
%
% OUTPUT: log_likelihood - log likelihood value

    N = length(hat_theta);
    
    % Initialize
    log_likelihood = 0;
    total_J = 0;
    
    % First part: constant term
    for i = 1:N
        J_i = length(hat_theta{i});
        total_J = total_J + J_i;
    end
    log_likelihood = log_likelihood - (N * total_J / 2) * log(2 * pi);
    
    % Second part: product of (nu + sigma^2_ij)^(-1/2)
    for i = 1:N
        J_i = length(hat_theta{i});
        for j = 1:J_i
            log_likelihood = log_likelihood - 0.5 * log(nu + rep_base_se2{i}(j));
        end
    end
    
    % Third part: exp term with sum_i sum_j
    sum_term = 0;
    for i = 1:N
        J_i = length(hat_theta{i});
        for j = 1:J_i
            sum_term = sum_term + (hat_theta{i}(j)^2) / (nu + rep_base_se2{i}(j));
        end
    end
    log_likelihood = log_likelihood - 0.5 * sum_term;
    
    % Fourth part: product of inverse sum terms
    for i = 1:N
        J_i = length(hat_theta{i});
        sum_inverse = 0;
        for j = 1:J_i
            sum_inverse = sum_inverse + 1 / (nu + rep_base_se2{i}(j));
        end
        log_likelihood = log_likelihood - 0.5 * log(sum_inverse);
    end
    
    % Fifth part: exp term with weighted sums
    for i = 1:N
        J_i = length(hat_theta{i});
        
        % Calculate the sum terms
        sum_inverse = 0;
        weighted_sum = 0;
        
        for j = 1:J_i
            sum_inverse = sum_inverse + 1 / (nu + rep_base_se2{i}(j));
            weighted_sum = weighted_sum + hat_theta{i}(j) / (nu + rep_base_se2{i}(j));
        end
        
        % Add to log likelihood
        log_likelihood = log_likelihood + 0.5 * (weighted_sum^2) / sum_inverse;
    end
end